Linear Regression
Nhóm 19
01-12-2024
Thư viện
library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errorslibrary(knitr)
library(janitor)##
## Attaching package: 'janitor'
##
## The following objects are masked from 'package:stats':
##
## chisq.test, fisher.testlibrary(stringr)
library(stats)
library(boot)
library(caret)## Loading required package: lattice
##
## Attaching package: 'lattice'
##
## The following object is masked from 'package:boot':
##
## melanoma
##
##
## Attaching package: 'caret'
##
## The following object is masked from 'package:purrr':
##
## liftlibrary(leaps)
library(glmnet)## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
##
## Loaded glmnet 4.1-8library(Metrics)##
## Attaching package: 'Metrics'
##
## The following objects are masked from 'package:caret':
##
## precision, recalllibrary(car)## Loading required package: carData
##
## Attaching package: 'car'
##
## The following object is masked from 'package:boot':
##
## logit
##
## The following object is masked from 'package:dplyr':
##
## recode
##
## The following object is masked from 'package:purrr':
##
## somelibrary(visreg)
library(corrplot)## corrplot 0.94 loadedlibrary(lmPerm)
library(dplyr)
library(mgcv)## Loading required package: nlme
##
## Attaching package: 'nlme'
##
## The following object is masked from 'package:dplyr':
##
## collapse
##
## This is mgcv 1.9-1. For overview type 'help("mgcv-package")'.Đọc và xử lí dữ liệu
data <- read.csv("cleaned_fifa_eda_stats.csv")
head(data)## id name age nationality overall potential
## 1 158023 L. Messi 31 Argentina 94 94
## 2 20801 Cristiano Ronaldo 33 Portugal 94 94
## 3 190871 Neymar Jr 26 Brazil 92 93
## 4 193080 De Gea 27 Spain 91 93
## 5 192985 K. De Bruyne 27 Belgium 91 92
## 6 183277 E. Hazard 27 Belgium 91 91
## club value wage preferred_foot international_reputation
## 1 FC Barcelona 110.5 0.565 Left 5
## 2 Juventus 77.0 0.405 Right 5
## 3 Paris Saint-Germain 118.5 0.290 Right 5
## 4 Manchester United 72.0 0.260 Right 4
## 5 Manchester City 102.0 0.355 Right 4
## 6 Chelsea 93.0 0.340 Right 4
## weak_foot skill_moves attacking_work_rate deffensive_work_rate body_type
## 1 4 4 Medium Medium Lean
## 2 4 5 High Low Stocky
## 3 5 5 High Medium Lean
## 4 3 1 Medium Medium Lean
## 5 5 4 High High Normal
## 6 4 4 High Medium Normal
## position jersey_number joined loaned_from contract_valid_until height
## 1 Forward 10 2004-07-01 Not on loan 2021-12-31 67
## 2 Forward 7 2018-07-10 Not on loan 2022-12-31 74
## 3 Forward 10 2017-08-03 Not on loan 2022-12-31 69
## 4 Goalkeeper 1 2011-07-01 Not on loan 2020-12-31 76
## 5 Midfielder 7 2015-08-30 Not on loan 2023-12-31 71
## 6 Forward 10 2012-07-01 Not on loan 2020-12-31 68
## weight crossing finishing heading_accuracy short_passing volleys dribbling
## 1 159 84 95 70 90 86 97
## 2 183 84 94 89 81 87 88
## 3 150 79 87 62 84 84 96
## 4 168 17 13 21 50 13 18
## 5 154 93 82 55 92 82 86
## 6 163 81 84 61 89 80 95
## curve fk_accuracy long_passing ball_control acceleration sprint_speed agility
## 1 93 94 87 96 91 86 91
## 2 81 76 77 94 89 91 87
## 3 88 87 78 95 94 90 96
## 4 21 19 51 42 57 58 60
## 5 85 83 91 91 78 76 79
## 6 83 79 83 94 94 88 95
## reactions balance shot_power jumping stamina strength long_shots aggression
## 1 95 95 85 68 72 59 94 48
## 2 96 70 95 95 88 79 93 63
## 3 94 84 80 61 81 49 82 56
## 4 90 43 31 67 43 64 12 38
## 5 91 77 91 63 90 75 91 76
## 6 90 94 82 56 83 66 80 54
## interceptions positioning vision penalties composure marking standing_tackle
## 1 22 94 94 75 96 33 28
## 2 29 95 82 85 95 28 31
## 3 36 89 87 81 94 27 24
## 4 30 12 68 40 68 15 21
## 5 61 87 94 79 88 68 58
## 6 41 87 89 86 91 34 27
## sliding_tackle gk_diving gk_handling gk_kicking gk_positioning gk_reflexes
## 1 26 6 11 15 14 8
## 2 23 7 11 15 14 11
## 3 33 9 9 15 15 11
## 4 13 90 85 87 88 94
## 5 51 15 13 5 10 13
## 6 22 11 12 6 8 8
## release_clause
## 1 226.5
## 2 127.1
## 3 228.1
## 4 138.6
## 5 196.4
## 6 172.1Ta loại bỏ các cột không đóng góp giá trị như id, name, nationality, loan_from, club
- id, name: do là các giá trị riêng biệt
- club : do các biến này có thể đóng góp trong thực tế (ví dụ như cầu thủ của Real Madrid thì tiềm năng hơn) tuy nhiên ta sẽ loại bỏ nó vì có quá nhiều câu lạc bộ sẽ làm tăng quá nhiều độ phức tạp của mô hình, tương tự với ‘nationality’ và ‘loan_from’.
data <- data|> dplyr::select(-c(id,name))
data <- data |> dplyr::select(-c(loaned_from,club,nationality))Xử lí dữ liệu:
- Ta chỉnh các biến định tính về dạng ‘factor’
- Các biến có dạng mm-dd-yyyy thì ta chỉ giữ lại năm (yyyy).
data_numeric <- data |> dplyr::select(where(is.numeric))
all_colname <- names(data)
category_col <- setdiff(all_colname, names(data_numeric))
for (i in 1:length(category_col)) {
data <- data |>
mutate(across(all_of(category_col[i]), as.factor))
}
data$contract_valid_until <- as.numeric(substr(data$contract_valid_until, 1, 4))
data$joined <- as.numeric(substr(data$joined, 1, 4))Tách dữ liệu làm hai tập train test
set.seed(100)
id_train <- sample(1:nrow(data),size = as.integer(0.8* nrow(data)))
id_test <- setdiff(1:nrow(data), id_train)
data_train <- data[id_train,]
data_test <- data[id_test,]I. Linear Regression cho Overall
1. Mô hình hồi quy với tất cả các biến
model <- lm(formula = overall ~ ., data = data_train)result_display <- function(model, data_test, target){
# Phần dư và các chỉ số trên tập huấn luyện
residuals_train <- residuals(model)
mae_train <- mean(abs(residuals_train))
mse_train <- mean(residuals_train^2)
adjr2_train <- summary(model)$adj.r.squared
# Dự đoán trên tập kiểm tra
y_predict <- predict(model, newdata = data_test)
# Tính các chỉ số trên tập kiểm tra
mse_test <- mean((data_test[[target]] - y_predict)^2)
mae_test <- mean(abs(data_test[[target]] - y_predict))
SSE <- sum((data_test[[target]] - y_predict)^2)
SST <- sum((data_test[[target]] - mean(data_test[[target]]))^2)
r_squared_test <- 1 - (SSE / SST)
n_test <- length(data_test[[target]]) # số mẫu trong tập kiểm tra
k <- length(model$coefficients) - 1 # số tham số trong mô hình (không tính hằng số)
adj_r_squared_test <- 1 - ((1 - r_squared_test) * (n_test - 1) / (n_test - k - 1))
# Tạo bảng kết quả
results <- data.frame(
Metric = c("MSE", "MAE", "Adjusted R2"),
Train = c(mse_train, mae_train, adjr2_train),
Test = c(mse_test, mae_test, adj_r_squared_test)
)
return(results)
}result_display(model,data_test,'overall')## Metric Train Test
## 1 MSE 3.105704 3.1407562
## 2 MAE 1.376685 1.3953978
## 3 Adjusted R2 0.936822 0.93349852. Lựa chọn mô hình
2.1 Step_wise với kFolds Cross-Validaiton
library(leaps)
set.seed(100)
k = 5
folds<-sample(rep(1:k,length=nrow(data_train))) predict.regsubsets <- function(object, newdata, id_model){
form <- as.formula(object$call[[2]])
x_mat <- model.matrix(form, newdata)
coef_est <- coef(object, id = id_model)
x_vars <- names(coef_est)
res <- x_mat[, x_vars] %*% coef_est
return(as.numeric(res))
}Ta sẽ dùng phương pháp là ‘backward’ để giảm thiểu thời gian cho việc chọn biến vì nếu như dùng ‘exhaustive’ cho quá nhiều biến thì rất tốn thời gian và thậm chí không tìm ra.
cv_error_fifa_rj<-matrix(0,nrow=k,ncol=dim(data_train)[2])
for (r in 1: k){
train <- data_train[folds !=r,]
validation <- data_train[folds ==r,]
out_subset_fifa_folds<-regsubsets(x=overall~ .,data=train,method="backward",nvmax=dim(train)[2])
for(j in 1:dim(train)[2]){
pred_rj<-predict(out_subset_fifa_folds, newdata=validation,id_model=j)
cv_error_fifa_rj[r,j]<-(mean((validation$overall-pred_rj)^2))
}
}
cv_error_fifa<-colMeans(cv_error_fifa_rj)ggplot(data= data.frame(x= c(1:dim(data_train)[2]),y=cv_error_fifa),aes(x=x,y=y)) +
geom_point() +
geom_line()+
labs(x="Number of predictors",y="RMSE") +
theme_bw()
cv_error_fifa## [1] 27.593196 7.339429 5.373930 5.149759 4.855027 4.404690 4.130525
## [8] 3.936668 3.884836 3.804003 3.717270 3.631348 3.550559 3.468331
## [15] 3.412661 3.349174 3.327005 3.318421 3.295094 3.278514 3.254145
## [22] 3.240575 3.229941 3.214501 3.208693 3.203589 3.198037 3.188483
## [29] 3.184845 3.174307 3.165740 3.158478 3.155494 3.149582 3.147407
## [36] 3.147207 3.145309 3.145252 3.143827 3.141378 3.139023 3.138328
## [43] 3.136577 3.136931 3.137426 3.137658 3.138497 3.139476 3.139936
## [50] 3.139705 3.139229 3.139325 3.139597Dưới đây là mô hình tốt nhất được lựa chọn từ phương pháp này
out_subset_fifa_2 <- regsubsets(x = overall ~ ., data = data_train,method = "backward",nvmax = dim(train)[2])
a <- names(coef(out_subset_fifa_2, id = which.min(cv_error_fifa)))data_train_subset <- data_train|> dplyr::select(-c(wage,preferred_foot,weak_foot,attacking_work_rate,deffensive_work_rate,body_type,position,joined,height,dribbling,curve,fk_accuracy,long_passing,jumping,long_shots,interceptions))Ta thấy các thông số đo lường tuy có giảm nhưng không quá đáng kể, nhưng ta đã giảm được độ phức tạp của thuật toán khá nhiều khi bỏ đi một số các biến đầu vào.
model <- lm(formula = overall ~ ., data = data_train_subset)
result_display(model,data_test,'overall')## Metric Train Test
## 1 MSE 3.2069351 3.2131129
## 2 MAE 1.3995615 1.4078213
## 3 Adjusted R2 0.9348659 0.93240052.2 Ridge và Lasso
x <- model.matrix(overall ~ ., data = data_train)[,-1]
y <- data_train$overalllibrary(glmnet)out_cv_lasso <- cv.glmnet(x = x, y = y, alpha = 1, type.measure = "mse", nfolds = 5,family = "gaussian")
lambda_grid <- 10^seq(from = 10, to =-2, length = 100)
beta_lambda_lasso <- out_cv_lasso$lambda.min
out_lasso_md <- glmnet(x = x, y = y, alpha = 1,lambda = lambda_grid, family = "gaussian")
predict(out_lasso_md, s = beta_lambda_lasso, type = "coefficients")## 58 x 1 sparse Matrix of class "dgCMatrix"
## s1
## (Intercept) -1.528883e+02
## age 4.889756e-01
## potential 4.681959e-01
## value 1.179531e-01
## wage .
## preferred_footRight .
## international_reputation -3.668166e-01
## weak_foot 5.964669e-03
## skill_moves 7.541796e-01
## attacking_work_rateLow 3.901313e-01
## attacking_work_rateMedium .
## deffensive_work_rate Low 3.059105e-01
## deffensive_work_rate Medium -3.518333e-02
## body_typeNormal 6.995978e-02
## body_typeStocky 9.750410e-02
## positionForward -1.001122e+00
## positionGoalkeeper .
## positionMidfielder -9.541762e-01
## jersey_number -9.447815e-03
## joined 5.150668e-04
## contract_valid_until 7.033457e-02
## height .
## weight 8.433614e-03
## crossing 2.241263e-03
## finishing 1.150634e-02
## heading_accuracy 3.790878e-02
## short_passing 4.478217e-02
## volleys .
## dribbling 4.584999e-05
## curve .
## fk_accuracy .
## long_passing 1.388231e-03
## ball_control 6.799523e-02
## acceleration 1.411227e-02
## sprint_speed 2.390074e-02
## agility 4.878108e-03
## reactions 1.347731e-01
## balance -1.506928e-02
## shot_power 1.439248e-02
## jumping 2.016065e-03
## stamina 2.921972e-02
## strength 2.308028e-02
## long_shots .
## aggression 1.696301e-03
## interceptions .
## positioning -1.391864e-02
## vision -5.738963e-03
## penalties -3.289357e-03
## composure 3.749191e-02
## marking 6.715413e-03
## standing_tackle 4.937880e-04
## sliding_tackle -6.323383e-05
## gk_diving 4.587257e-02
## gk_handling 3.642575e-02
## gk_kicking 7.611226e-03
## gk_positioning 2.387573e-02
## gk_reflexes 5.079425e-02
## release_clause .Cross-Validation cho mô hình đã bỏ đi các biến theo Lasso
data_train_lasso <- data_train|>dplyr::select(-c(release_clause,interceptions,long_shots,curve,fk_accuracy,volleys,weight,attacking_work_rate,wage,preferred_foot))
mse_list <- c()
# kfolds = 5
for (i in 1: 5){
train <- data_train_lasso[folds != i,]
validate <- data_train_lasso[folds == i,]
fifa_lasso <- lm(formula = overall ~.,data = train)
mse <- (mean((validate$overall-predict(newdata = validate,fifa_lasso))^2))
mse_list <- c(mse_list,mse)
}mse_cv <- mean(mse_list)
cat("RMSE: ",mse_cv)## RMSE: 3.160872Các thông số bằng phương pháp này cho kết quả tốt hơn so với phương pháp trước nên vì thế ta sẽ dùng mô hình này cho các bước tiếp theo.
data_train_lasso <- data_train|>dplyr::select(-c(release_clause,interceptions,long_shots,curve,fk_accuracy,volleys,weight,attacking_work_rate,wage,preferred_foot))
model <-lm(formula = overall ~.,data = data_train_lasso)
result_display(model,data_test,'overall')## Metric Train Test
## 1 MSE 3.1328745 3.178213
## 2 MAE 1.3833992 1.402746
## 3 Adjusted R2 0.9363221 0.9329313. Chuẩn đoán mô hình
3.1 Kiểm tra tính tuyến tính của mô hình
model <- lm(data = data_train_lasso, formula = overall ~. )library(ggplot2)
ggplot(model,aes(x = .fitted, y = .resid))+
geom_point()+
geom_smooth(method = 'loess',se = F)+
geom_hline(yintercept = 0,linetype = 'dashed')+
theme_bw()+
labs(x = "Fitted values", y = "Residuals")## `geom_smooth()` using formula = 'y ~ x'
Ta thấy dựa vào hình thì giá trị thặng dư khá tốt ở lúc đầu tuy nhiên phần cuối lại có dấu hiệu là đường cong nên vì thế ta sẽ hiệu chỉnh sau.
3.2 Kiểm tra tính tuyến tính từng phần
terms_md <- predict(model, type = "terms")
partial_resid_md <- residuals(model, type = "partial")
col = names(data.frame(terms_md))
# Duyệt qua tất cả các cột trong data và vẽ biểu đồ cho từng cột
for (i in 1:length(col)){
# Tạo dataframe cho mỗi cột
data_part_resid_md <- tibble(
data_ = data_train[, col[i]],
terms = terms_md[, col[i]],
partial_resid = partial_resid_md[, col[i]]
)
# Vẽ biểu đồ và hiển thị biểu đồ cho từng cột
p <- ggplot(data_part_resid_md, mapping = aes(data_, partial_resid)) +
geom_point() +
geom_smooth(method = "loess", se = FALSE, linetype = "dashed", color = "forestgreen") +
geom_line(aes(x = data_, y = terms), color = "blue") +
labs(x = col[i], y = "Partial Residuals") +
theme_bw()
# In ra biểu đồ
print(p)
}## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : at 0.98## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : radius 0.0004## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : all data on boundary of neighborhood. make span bigger## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 0.98## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 0.02## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 1## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : zero-width neighborhood. make span bigger## Warning: Failed to fit group -1.
## Caused by error in `predLoess()`:
## ! NA/NaN/Inf in foreign function call (arg 5)
## `geom_smooth()` using formula = 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 0.98## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 2.02## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 5.9127e-15## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : There are other near singularities as well. 4.0804
## `geom_smooth()` using formula = 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 2## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 1## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 0## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : There are other near singularities as well. 1
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 2017## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 1## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 0
## `geom_smooth()` using formula = 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 2020## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 1## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 0
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
3.3 Kiểm tra đồng nhất phương sai
ggplot(model, aes(.fitted, sqrt(abs(.stdresid)))) +
geom_point(na.rm = TRUE) +
geom_smooth(method = "loess", na.rm = TRUE, se = FALSE) +
labs(x = "Fitted Values", y = expression(sqrt("|Standardized residuals|"))) +
theme_bw()## `geom_smooth()` using formula = 'y ~ x'
Từ hình ảnh ta thấy các điểm đầu khá chuẩn chỉ có những điểm dữ lệu cuối có xu hướng đường cong rõ rệt tương tự như phần 3.1
3.4 Kiểm tra điểm ngoại lai
ggplot(model, aes(.hat, .stdresid)) +
geom_point(aes(size = .cooksd)) +
xlab("Leverage") + ylab("Standardized Residuals") +
scale_size_continuous("Cook's Distance", range = c(1, 6)) +
theme_bw() +
theme(legend.position = "bottom")
std_resid_md <- rstandard(model)
hat_values_md <- hatvalues(model)
cooks_D_md <- cooks.distance(model) data_cooks_md <- tibble(id_point = 1:nrow(data_train_lasso),
rstand = std_resid_md, hats = hat_values_md,
cooks = cooks_D_md, overall = data_train_lasso$overall)
data_cooks_md <- data_cooks_md |> arrange(desc(cooks))
head(data_cooks_md)## # A tibble: 6 × 5
## id_point rstand hats cooks overall
## <int> <dbl> <dbl> <dbl> <int>
## 1 879 -4.38 0.0399 0.0170 94
## 2 6110 -4.47 0.0336 0.0148 91
## 3 8803 -5.28 0.0239 0.0145 91
## 4 4357 -4.89 0.0195 0.0101 94
## 5 13229 -3.52 0.0224 0.00606 91
## 6 3209 -4.81 0.0106 0.00528 56Ta sẽ thực hiện loại bỏ Outlier cho mô hình
Những điểm nào có cook distance lớn hơn 4/n với n là tổng số quan sát thì ta sẽ loại bỏ những điểm này
id_exclude <- (data_cooks_md |> filter(cooks > 4/nrow(data_train_lasso)))$id_point
id_rest <- setdiff(data_cooks_md$id_point,id_exclude)
data_train_lasso <- data_train_lasso[id_rest,]model <- lm(data = data_train_lasso, formula = overall ~. )
result_display(model,data_test,'overall')## Metric Train Test
## 1 MSE 2.2461048 3.2155173
## 2 MAE 1.2121397 1.3954416
## 3 Adjusted R2 0.9514614 0.9321437Ta thấy mô hình cho kết quả tốt hơn ở tập train so với ban đầu. Vì các điểm Outlier đã bị loại bỏ giúp cho mô hình học tốt hơn. Tuy nhiên ở tập test đã tăng lên. Ta sẽ giải quyết vấn đề này ở phần 3
3.5 Kiểm tra đa cộng tuyến
model <- lm(data = data_train_lasso, formula = overall ~. )
library(car)
df_vif <- data.frame(vif(model))|>arrange(desc(GVIF)) |>filter(GVIF >5)
df_vif## GVIF Df GVIF..1..2.Df..
## position 340.117025 3 2.642032
## gk_diving 29.814994 1 5.460311
## gk_reflexes 29.737673 1 5.453226
## gk_handling 27.429321 1 5.237301
## gk_positioning 26.026394 1 5.101607
## standing_tackle 24.929892 1 4.992984
## gk_kicking 23.005440 1 4.796399
## sliding_tackle 22.969645 1 4.792666
## ball_control 19.398961 1 4.404425
## dribbling 16.015987 1 4.001998
## short_passing 13.112069 1 3.621059
## finishing 9.762605 1 3.124517
## positioning 9.463701 1 3.076313
## acceleration 8.902934 1 2.983779
## long_passing 7.514671 1 2.741290
## sprint_speed 7.368550 1 2.714507
## marking 6.802705 1 2.608200
## heading_accuracy 6.732468 1 2.594700
## crossing 6.111110 1 2.472066
## shot_power 5.123720 1 2.263564remove_multicollinearity <- function(feature, data_train) {
return (data_train |> dplyr::select(-all_of(feature)))
}while(TRUE){
model <- lm(data = data_train_lasso, formula = overall ~. )
df_vif <- data.frame(vif(model))|>arrange(desc(GVIF)) |>filter(GVIF >5)|>slice_max(GVIF, n = 1)
vec_feature <- row.names(df_vif)
if (length(vec_feature) == 0){
break
}
data_train_lasso <- remove_multicollinearity(vec_feature,data_train_lasso)
}model <- lm(data = data_train_lasso, formula = overall ~. )
result_display(model,data_test,'overall')## Metric Train Test
## 1 MSE 3.0568351 4.060881
## 2 MAE 1.3910795 1.550162
## 3 Adjusted R2 0.9340255 0.914720Sau khi loại bỏ hết đa cộng tuyến thì mô hình trở nên rất tê, tuy nhiên nó đã giảm được độ phức tạp cuẩ thuật toán khá nhiều
4. Mở rộng mô hình
upgrade_model <- lm(formula = overall ~ poly(age,4)+poly(value,2) +poly(potential,2)+., data = data_train_lasso)
result_display(upgrade_model,data_test,'overall')## Metric Train Test
## 1 MSE 1.5574229 2.6839806
## 2 MAE 0.9411552 1.0527663
## 3 Adjusted R2 0.9663734 0.9434984Ta sẽ sử dụng hàm mũ cho các biến như ‘age’, ‘value’, ‘potential’ ta thấy mô hình đã cải thiện rõ rệt và cho thông số rất tốt.
data_tst <- data_test|>dplyr::select(names(data_train_lasso))
new_obs <- data_tst[10,]
new_obs## age overall potential value international_reputation weak_foot skill_moves
## 57 26 86 87 52 3 4 4
## deffensive_work_rate body_type jersey_number joined contract_valid_until
## 57 Medium Lean 10 2016 2023
## height crossing heading_accuracy long_passing sprint_speed agility reactions
## 57 69 73 62 71 93 91 86
## balance shot_power jumping stamina strength aggression vision penalties
## 57 86 82 75 84 67 73 82 71
## composure marking
## 57 80 425. Xây dựng khoảng tin cậy
5.1 Xay dựng khoảng tin cậy cho từng biến
fun_boot_md <- function(data, ind, formula, ...){
data_new <- data[ind,]
out_md <- lm(formula = formula, data = data_new, ...)
return(out_md$coefficients)
}
set.seed(100)
out_boot_md <- boot(data = data_train_lasso, statistic = fun_boot_md, R = 1000,
formula = overall ~ poly(age,4)+ poly(value,2)+ poly(potential,2)+ .,parallel = "multicore")pvals <- sapply(1:ncol(out_boot_md$t),function(x) {
qt0 <- mean(out_boot_md$t[, x] <= 0)
if (!is.na(qt0)){
if (qt0 < 0.5) {
return(2*qt0)
} else {
return(2*(1- qt0))
}
}
}
)
pvals <- unlist(pvals)
cat(pvals, sep = " ")## 0.674 0 0 0 0 0 0 0 0 0 0.386 0 0.102 0.538 0 0 0 0.23 0 0.002 0 0 0 0 0 0 0 0.072 0 0.748 0 0.164 0.112 0 0 0Bên trên và p_value cho từng biến
5.2 Khoảng tin cậy cho trung bình
fun_model_predict_mu <- function(data,ind,formula,new_observation){
data_new <- data[ind,]
model <- lm(formula = formula,data = data_new)
pred <- predict(model,newdata = new_observation)
return (pred)
}library(boot)
set.seed(100)
out_boot_pred <- boot(data = data_train_lasso, statistic = fun_model_predict_mu,formula = overall ~ poly(age,4)+ poly(value,2)+ poly(potential,2)+.,new_observation = new_obs, R = 1000, parallel = "multicore")df <- data.frame(t = out_boot_pred$t)
# Create the histogram using ggplot
ggplot(df, aes(x = t)) +
geom_histogram(binwidth = 0.2, fill = "lightblue", color = "black") +
labs(x = "Bootstrap Estimates", y = "Frequency") +
theme_minimal()
quantile(out_boot_pred$t,c(0.025,0.975))## 2.5% 97.5%
## 82.71068 85.240965.3 Khoảng tin cậy cho biến dự đoán
resid <- residuals(upgrade_model)
y_pd_pci <- out_boot_pred$t + sample(resid, size = 1000, replace = TRUE)
quantile(y_pd_pci, probs = c(0.025, 0.975))## 2.5% 97.5%
## 81.41045 86.92959II. Linear Regression cho Potential
1. Mô hình hồi quy với tất cả các biến
model <- lm(formula = potential ~ ., data = data_train)result_display(model,data_test,'potential')## Metric Train Test
## 1 MSE 5.9263871 5.8796982
## 2 MAE 1.9126669 1.9026382
## 3 Adjusted R2 0.8443281 0.83492532. Lựa chọn mô hình
2.1 Step_wise và kFolds Cross-Validaiton
library(leaps)
set.seed(101)
k = 5
folds<-sample(rep(1:k,length=nrow(data_train))) predict.regsubsets <- function(object, newdata, id_model){
form <- as.formula(object$call[[2]])
x_mat <- model.matrix(form, newdata)
coef_est <- coef(object, id = id_model)
x_vars <- names(coef_est)
res <- x_mat[, x_vars] %*% coef_est
return(as.numeric(res))
}cv_error_fifa_rj<-matrix(0,nrow=k,ncol=dim(data_train)[2])
for (r in 1: k){
train <- data_train[folds !=r,]
validation <- data_train[folds ==r,]
out_subset_fifa_folds<-regsubsets(x=potential~ .,data=train,method="backward",nvmax=dim(train)[2])
for(j in 1:dim(train)[2]){
pred_rj<-predict(out_subset_fifa_folds, newdata=validation,id_model=j)
cv_error_fifa_rj[r,j]<-(mean((validation$potential-pred_rj)^2))
}
}
cv_error_fifa<-colMeans(cv_error_fifa_rj)ggplot(data= data.frame(x= c(1:dim(data_train)[2]),y=cv_error_fifa),aes(x=x,y=y)) +
geom_point() +
geom_line()+
labs(x="Number of predictors",y="RMSE") +
theme_bw()
out_subset_fifa_2 <- regsubsets(x = potential ~ ., data = data_train,method = "backward",nvmax = dim(data_train)[2])
a <- names(coef(out_subset_fifa_2, id = which.min(cv_error_fifa)))setdiff(names(data_train),a)## [1] "potential" "wage" "preferred_foot"
## [4] "skill_moves" "attacking_work_rate" "deffensive_work_rate"
## [7] "body_type" "position" "weight"
## [10] "heading_accuracy" "short_passing" "volleys"
## [13] "dribbling" "curve" "fk_accuracy"
## [16] "ball_control" "reactions" "shot_power"
## [19] "jumping" "aggression" "interceptions"
## [22] "gk_diving" "gk_kicking"data_train_subset <- data_train|> dplyr::select(-c(
wage, preferred_foot, skill_moves, attacking_work_rate,
deffensive_work_rate, body_type, position, weight, heading_accuracy,
dribbling, curve, fk_accuracy, reactions, shot_power, jumping,
interceptions, gk_kicking)
)model <- lm(formula = potential ~ ., data = data_train_subset)
result_display(model,data_test,'potential')## Metric Train Test
## 1 MSE 5.9605624 5.9274720
## 2 MAE 1.9178459 1.9110744
## 3 Adjusted R2 0.8436898 0.83469582.2 Ridge và Lasso
x <- model.matrix(potential ~ ., data = data_train)[,-1]
y <- data_train$potentiallibrary(glmnet)out_cv_lasso <- cv.glmnet(x = x, y = y, alpha = 1, type.measure = "mse", nfolds = 5,family = "gaussian")
lambda_grid <- 10^seq(from = 10, to =-2, length = 100)
beta_lambda_lasso <- out_cv_lasso$lambda.min
out_lasso_md <- glmnet(x = x, y = y, alpha = 1,lambda = lambda_grid, family = "gaussian")
predict(out_lasso_md, s = beta_lambda_lasso, type = "coefficients")## 58 x 1 sparse Matrix of class "dgCMatrix"
## s1
## (Intercept) -3.638345e+01
## age -9.240860e-01
## overall 8.707526e-01
## value -1.681563e-02
## wage .
## preferred_footRight -3.750818e-02
## international_reputation 1.002604e+00
## weak_foot 3.983963e-02
## skill_moves .
## attacking_work_rateLow 1.120708e-02
## attacking_work_rateMedium -8.281785e-04
## deffensive_work_rate Low -1.386661e-01
## deffensive_work_rate Medium -6.593011e-02
## body_typeNormal -2.084916e-01
## body_typeStocky -2.174540e-01
## positionForward 3.466959e-01
## positionGoalkeeper -5.170995e-01
## positionMidfielder -1.499088e-02
## jersey_number 1.200514e-02
## joined -3.078260e-02
## contract_valid_until 6.910695e-02
## height -3.453569e-02
## weight -2.620536e-03
## crossing -1.238857e-02
## finishing -1.430527e-04
## heading_accuracy .
## short_passing .
## volleys -1.254810e-03
## dribbling .
## curve .
## fk_accuracy -2.195305e-04
## long_passing -2.264711e-03
## ball_control .
## acceleration -4.019492e-03
## sprint_speed -1.835308e-02
## agility -3.280621e-03
## reactions .
## balance 7.270086e-03
## shot_power .
## jumping .
## stamina -4.183314e-02
## strength -1.330754e-02
## long_shots -7.347920e-03
## aggression .
## interceptions .
## positioning -6.392837e-03
## vision 1.304703e-02
## penalties 1.525798e-02
## composure 2.325342e-02
## marking 8.231183e-03
## standing_tackle .
## sliding_tackle 7.824327e-03
## gk_diving .
## gk_handling .
## gk_kicking .
## gk_positioning 3.407980e-03
## gk_reflexes -3.800857e-03
## release_clause .Cross- Validation cho Lasso Model
data_train_lasso <- data_train|>dplyr::select(-c(release_clause,gk_reflexes,gk_kicking,gk_handling,gk_diving,sliding_tackle,interceptions,shot_power,reactions,ball_control,long_passing,fk_accuracy,curve,dribbling,short_passing,finishing,skill_moves,wage))
mse_list <- c()
# kfolds = 5
for (i in 1: 5){
train <- data_train_lasso[folds != i,]
validate <- data_train_lasso[folds == i,]
fifa_lasso <- lm(formula = potential ~.,data = train)
mse <- (mean((validate$potential-predict(newdata = validate,fifa_lasso))^2))
mse_list <- c(mse_list,mse)
}mse_cv <- mean(mse_list)
cat("RMSE: ",mse_cv)## RMSE: 6.015211model <-lm(formula = potential ~.,data = data_train_lasso)
result_display(model,data_test,'potential')## Metric Train Test
## 1 MSE 5.9694848 5.9365298
## 2 MAE 1.9179766 1.9108106
## 3 Adjusted R2 0.8434086 0.8342419Ta chọn bộ data ‘data_train_subset’ cho những phần tiếp theo.
3. Chuẩn đoán mô hình
3.1 Kiểm tra tính tuyến tính của mô hình
model <- lm(data = data_train_subset, formula = potential ~. )library(ggplot2)
ggplot(model,aes(x = .fitted, y = .resid))+
geom_point()+
geom_smooth(method = 'loess',se = F)+
geom_hline(yintercept = 0,linetype = 'dashed')+
theme_bw()+
labs(x = "Fitted values", y = "Residuals")## `geom_smooth()` using formula = 'y ~ x'
Tính tuyến tính của mô hình là không thỏa chủ yếu ở những điểm đầu và cuối.
3.2 Kiểm tra tính tuyến tính từng phần
terms_md <- predict(model, type = "terms")
partial_resid_md <- residuals(model, type = "partial")
col = names(data.frame(terms_md))
# Duyệt qua tất cả các cột trong data và vẽ biểu đồ cho từng cột
for (i in 1:length(col)){
# Tạo dataframe cho mỗi cột
data_part_resid_md <- tibble(
data_ = data_train[, col[i]],
terms = terms_md[, col[i]],
partial_resid = partial_resid_md[, col[i]]
)
# Vẽ biểu đồ và hiển thị biểu đồ cho từng cột
p <- ggplot(data_part_resid_md, mapping = aes(data_, partial_resid)) +
geom_point() +
geom_smooth(method = "loess", se = FALSE, linetype = "dashed", color = "forestgreen") +
geom_line(aes(x = data_, y = terms), color = "blue") +
labs(x = col[i], y = "Partial Residuals") +
theme_bw()
# In ra biểu đồ
print(p)
}## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : at 0.98## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : radius 0.0004## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : all data on boundary of neighborhood. make span bigger## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 0.98## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 0.02## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 1## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : zero-width neighborhood. make span bigger## Warning: Failed to fit group -1.
## Caused by error in `predLoess()`:
## ! NA/NaN/Inf in foreign function call (arg 5)
## `geom_smooth()` using formula = 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 0.98## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 2.02## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 5.9127e-15## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : There are other near singularities as well. 4.0804
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 2017## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 1## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 0
## `geom_smooth()` using formula = 'y ~ x'## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : pseudoinverse used at 2020## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : neighborhood radius 1## Warning in simpleLoess(y, x, w, span, degree = degree, parametric = parametric,
## : reciprocal condition number 0
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
3.3 Kiểm tra đồng nhất phương sai
ggplot(model, aes(.fitted, sqrt(abs(.stdresid)))) +
geom_point(na.rm = TRUE) +
geom_smooth(method = "loess", na.rm = TRUE, se = FALSE) +
labs(x = "Fitted Values", y = expression(sqrt("|Standardized residuals|"))) +
theme_bw()## `geom_smooth()` using formula = 'y ~ x'
Mô hình không đồng nhất phương sai
3.4 Kiểm tra điểm ngoại lai
ggplot(model, aes(.hat, .stdresid)) +
geom_point(aes(size = .cooksd)) +
xlab("Leverage") + ylab("Standardized Residuals") +
scale_size_continuous("Cook's Distance", range = c(1, 6)) +
theme_bw() +
theme(legend.position = "bottom")
std_resid_md <- rstandard(model)
hat_values_md <- hatvalues(model)
cooks_D_md <- cooks.distance(model)Ta sẽ thực hiện loại bỏ Outlier cho mô hình
Những điểm nào có cook distance lớn hơn 4/n với n là tổng số quan sát thì ta sẽ loại bỏ những điểm này
data_cooks_md <- tibble(id_point = 1:nrow(data_train_subset),
rstand = std_resid_md, hats = hat_values_md,
cooks = cooks_D_md)
data_cooks_md <- data_cooks_md |> arrange(desc(cooks))
head(data_cooks_md)## # A tibble: 6 × 4
## id_point rstand hats cooks
## <int> <dbl> <dbl> <dbl>
## 1 4357 3.41 0.0434 0.0147
## 2 8052 2.22 0.0503 0.00724
## 3 12025 1.79 0.0461 0.00429
## 4 3485 3.47 0.0105 0.00355
## 5 2317 5.17 0.00284 0.00211
## 6 1716 4.49 0.00368 0.00206id_exclude <- (data_cooks_md |> filter(cooks > 4/nrow(data_train_subset)))$id_pointid_rest <- setdiff(data_cooks_md$id_point,id_exclude)data_train_subset <- data_train_subset[id_rest,]model <- lm(data = data_train_subset, formula = potential ~. )
result_display(model,data_test,'potential')## Metric Train Test
## 1 MSE 4.3800636 6.0144317
## 2 MAE 1.6926823 1.9039860
## 3 Adjusted R2 0.8786617 0.8322707Ta thấy mô hình cho kết quả tốt hơn ở tập train so với ban đầu. Vì các điểm Outlier đã bị loại bỏ giúp cho mô hình học tốt hơn. Tuy nhiên ở tập test đã tăng lên. Ta sẽ giải quyết vấn đề này ở phần 3
3.5 Kiểm tra đa cộng tuyến
model <- lm(data = data_train_subset, formula = potential ~. )
library(car)
data.frame(vif = vif(model))|>arrange(desc(vif))## vif
## value 101.815814
## release_clause 95.415415
## gk_reflexes 27.058978
## gk_diving 26.881912
## gk_positioning 24.627297
## gk_handling 24.619208
## standing_tackle 24.408502
## sliding_tackle 22.429390
## ball_control 16.145992
## short_passing 13.199222
## finishing 10.640687
## positioning 8.979948
## acceleration 8.723754
## long_passing 7.297106
## sprint_speed 7.285258
## long_shots 7.202759
## marking 6.633039
## overall 6.346339
## volleys 6.256924
## crossing 5.177654
## agility 4.672618
## balance 4.498720
## penalties 4.434403
## vision 4.170814
## stamina 3.830155
## aggression 3.630437
## composure 3.620232
## height 3.391360
## strength 2.776957
## international_reputation 2.088940
## age 2.011132
## weak_foot 1.178769
## contract_valid_until 1.173539
## jersey_number 1.108738
## joined 1.103443while(TRUE){
model <- lm(data = data_train_subset, formula = potential ~. )
df_vif <- data.frame(vif = vif(model))|>arrange(desc(vif))|>filter(vif >5)|>slice_max(vif, n = 1)
vec_feature <- row.names(df_vif)
if (length(vec_feature) == 0){
break
}
data_train_subset<- remove_multicollinearity(vec_feature,data_train_subset)
}model <- lm(data = data_train_subset, formula = potential ~. )
result_display(model,data_test,'potential')## Metric Train Test
## 1 MSE 4.450271 6.1205107
## 2 MAE 1.705565 1.9174774
## 3 Adjusted R2 0.876844 0.82998364. Mở rộng mô hình
upgrade_model <- lm(formula = potential ~ poly(age,4)+., data = data_train_subset)
result_display(upgrade_model,data_test,'potential')## Metric Train Test
## 1 MSE 2.1285713 3.2472175
## 2 MAE 1.0426589 1.1800429
## 3 Adjusted R2 0.9410803 0.9096891summary <- summary(upgrade_model)
cat("Adjusted R2 train:",summary$adj.r.squared)## Adjusted R2 train: 0.9410803data_tst <- data_test|>dplyr::select(names(data_train_subset))
new_obs <- data_tst[4,]
new_obs## age overall potential international_reputation weak_foot jersey_number
## 39 33 88 88 3 2 1
## joined contract_valid_until height crossing volleys long_passing
## 39 2012 2021 76 12 12 34
## sprint_speed agility balance stamina strength aggression vision penalties
## 39 55 47 36 41 71 25 41 23
## composure marking release_clause
## 39 69 25 515. Xây dựng khoảng tin cậy
5.1 Xay dựng khoảng tin cậy cho từng biến
fun_boot_md <- function(data, ind, formula, ...){
data_new <- data[ind,]
out_md <- lm(formula = formula, data = data_new, ...)
return(out_md$coefficients)
}
set.seed(100)
out_boot_md <- boot(data = data_train_subset, statistic = fun_boot_md, R = 1000,
formula = potential ~ poly(age,4)+ .,parallel = "multicore")pvals <- sapply(1:ncol(out_boot_md$t),function(x) {
qt0 <- mean(out_boot_md$t[, x] <= 0)
if (!is.na(qt0)){
if (qt0 < 0.5) {
return(2*qt0)
} else {
return(2*(1- qt0))
}
}
}
)
pvals <- unlist(pvals)
cat(pvals, sep = " ")## 0.09 0 0 0 0 0 0 0.508 0 0 0 0.898 0 0.006 0.982 0.014 0.022 0.028 0 0.2 0 0.47 0.002 0.35 0 0.3685.2 Khoảng tin cậy cho trung bình
fun_model_predict_mu <- function(data,ind,formula,new_observation,knot){
data_new <- data[ind,]
model <- lm(formula = formula,data = data_new)
pred <- predict(model,newdata = new_observation)
return (pred)
}library(boot)
set.seed(100)
out_boot_pred <- boot(data = data_train_subset, statistic = fun_model_predict_mu,formula =
potential ~ poly(age,4)+.,new_observation = new_obs, R = 1000, parallel = "multicore")df <- data.frame(t = out_boot_pred$t)
# Create the histogram using ggplot
ggplot(df, aes(x = t)) +
geom_histogram(binwidth = 0.02, fill = "lightblue", color = "black") +
labs(x = "Bootstrap Estimates", y = "Frequency") +
theme_minimal()
quantile(out_boot_pred$t,c(0.025,0.975))## 2.5% 97.5%
## 88.19318 88.501225.3 Khoảng tin cậy cho biến dự đoán
resid <- residuals(upgrade_model)
y_pd_pci <- out_boot_pred$t + sample(resid, size = 1000, replace = TRUE)
quantile(y_pd_pci, probs = c(0.025, 0.975))## 2.5% 97.5%
## 85.18643 92.16681III. Linear Regresion cho Value
simplified_data <- data|> dplyr::mutate(wage = wage * 1000)value_lm <- lm(value ~ ., data = simplified_data)
summary(value_lm)##
## Call:
## lm(formula = value ~ ., data = simplified_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.5958 -0.1200 -0.0012 0.0994 12.2424
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.212e+00 8.374e+00 0.742 0.458256
## age -1.910e-02 2.186e-03 -8.737 < 2e-16 ***
## overall 3.463e-02 2.492e-03 13.894 < 2e-16 ***
## potential -1.887e-02 1.812e-03 -10.415 < 2e-16 ***
## wage 6.375e-03 4.121e-04 15.471 < 2e-16 ***
## preferred_footRight -2.878e-03 1.115e-02 -0.258 0.796337
## international_reputation 3.996e-01 1.621e-02 24.657 < 2e-16 ***
## weak_foot 1.405e-02 7.369e-03 1.906 0.056669 .
## skill_moves 7.028e-03 1.206e-02 0.583 0.560214
## attacking_work_rateLow -2.762e-04 2.425e-02 -0.011 0.990913
## attacking_work_rateMedium -8.047e-03 1.164e-02 -0.691 0.489405
## deffensive_work_rate Low -2.003e-03 2.107e-02 -0.095 0.924250
## deffensive_work_rate Medium 1.367e-03 1.309e-02 0.104 0.916853
## body_typeNormal 1.765e-02 1.038e-02 1.700 0.089120 .
## body_typeStocky 3.999e-02 2.122e-02 1.884 0.059518 .
## positionForward 2.842e-02 2.551e-02 1.114 0.265218
## positionGoalkeeper 4.991e-01 9.608e-02 5.194 2.08e-07 ***
## positionMidfielder 1.527e-02 1.795e-02 0.851 0.394799
## jersey_number -5.115e-04 2.921e-04 -1.751 0.079930 .
## joined 1.660e-03 2.223e-03 0.746 0.455377
## contract_valid_until -4.942e-03 3.711e-03 -1.332 0.183021
## height -4.716e-03 3.468e-03 -1.360 0.173931
## weight 7.156e-04 5.250e-04 1.363 0.172832
## crossing -2.233e-03 6.336e-04 -3.525 0.000425 ***
## finishing -3.162e-03 7.689e-04 -4.112 3.94e-05 ***
## heading_accuracy -9.633e-04 6.773e-04 -1.422 0.154999
## short_passing -2.754e-03 1.085e-03 -2.539 0.011135 *
## volleys 3.783e-03 6.564e-04 5.763 8.39e-09 ***
## dribbling -6.404e-04 9.505e-04 -0.674 0.500492
## curve -1.662e-03 6.483e-04 -2.564 0.010364 *
## fk_accuracy 1.958e-03 5.756e-04 3.402 0.000672 ***
## long_passing 2.529e-03 7.985e-04 3.168 0.001540 **
## ball_control -1.806e-03 1.184e-03 -1.525 0.127287
## acceleration 6.035e-04 8.892e-04 0.679 0.497316
## sprint_speed -9.413e-05 8.250e-04 -0.114 0.909165
## agility 7.978e-04 6.603e-04 1.208 0.226968
## reactions 1.772e-03 1.006e-03 1.761 0.078335 .
## balance -6.487e-04 6.748e-04 -0.961 0.336393
## shot_power -8.242e-04 6.693e-04 -1.231 0.218235
## jumping -2.772e-04 4.665e-04 -0.594 0.552348
## stamina 4.075e-03 5.562e-04 7.326 2.47e-13 ***
## strength 1.597e-04 6.561e-04 0.243 0.807728
## long_shots -7.718e-04 7.095e-04 -1.088 0.276655
## aggression -9.097e-04 5.067e-04 -1.795 0.072631 .
## interceptions 7.142e-04 7.226e-04 0.988 0.322936
## positioning 5.412e-04 7.098e-04 0.763 0.445760
## vision 1.383e-03 6.627e-04 2.087 0.036870 *
## penalties 1.828e-05 6.238e-04 0.029 0.976629
## composure -3.584e-05 7.386e-04 -0.049 0.961303
## marking 1.281e-03 5.884e-04 2.177 0.029521 *
## standing_tackle -2.176e-03 1.077e-03 -2.021 0.043266 *
## sliding_tackle -1.745e-03 1.003e-03 -1.739 0.082007 .
## gk_diving -2.558e-03 1.385e-03 -1.847 0.064759 .
## gk_handling -4.740e-03 1.394e-03 -3.400 0.000675 ***
## gk_kicking -1.530e-03 1.291e-03 -1.185 0.236041
## gk_positioning -1.812e-03 1.346e-03 -1.346 0.178178
## gk_reflexes -1.730e-03 1.369e-03 -1.264 0.206102
## release_clause 4.849e-01 8.764e-04 553.305 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5702 on 16585 degrees of freedom
## Multiple R-squared: 0.9901, Adjusted R-squared: 0.9901
## F-statistic: 2.909e+04 on 57 and 16585 DF, p-value: < 2.2e-16Nhận xét
- Dựa trên kết quả từ hàm summary, ta nhận thấy được độ chính xác của mô hình rất cao R-squared: 0.9901, Adjusted R-squared: 0.9901.
- Với mức ý nghĩa alpha = 0.0001 và kết quả từ cột Pr(>|t|), các biến giải thích có ảnh hưởng đáng kể đến mô hình bao gồm: age, overall, potential, wage, international_reputation, positionalGoalKeeper, crossing, finishing, volleys, fk_accuracy, stamina, gk_kicking, release_clause. Tuy nhiên, cần lưu ý rằng các giá trị p-value này dựa trên các giả định của mô hình hồi quy tuyến tính.
- Để đảm bảo tính chính xác và độ tin cậy của các hệ số ước lượng, nên áp dụng phương pháp bootstrap để kiểm định các hệ số thay vì chỉ dựa vào giả định ban đầu của mô hình.
1. Sự suy luận cho mô hình
- Vì dữ liệu có kích thước rất lớn, thời gian thực hiện Boostrap sẽ
rất lâu, nên ta thiết lập tham số
parallel = multicoređể thực hiện song song.
# Thực hiện phương pháp boostrap
fun_boot_md <- function(data, ind, formula, ...){
data_new <- data[ind,]
out_md <- lm(formula = formula, data = data_new, ...)
return(setNames(out_md$coefficients, names(out_md$coefficients)))
}
set.seed(84)
boot_value_lm <- boot(data = simplified_data, statistic = fun_boot_md, R = 1000, formula = value ~ ., parallel = "multicore")1.1 Kiểm định hệ số bằng phương pháp Bootstrap
pvals_boot_value <- sapply(1:ncol(boot_value_lm$t),function(x) {
qt0 <- mean(boot_value_lm$t[, x] <= 0)
if (qt0 < 0.5) {
return(2*qt0)
} else {
return(2*(1- qt0))
}
})
names(pvals_boot_value) <- names(boot_value_lm$t0)
pvals_boot_value[pvals_boot_value < 0.0001]## age overall potential
## 0 0 0
## wage international_reputation positionGoalkeeper
## 0 0 0
## crossing finishing volleys
## 0 0 0
## long_passing stamina gk_handling
## 0 0 0
## release_clause
## 0Nhận xét
- Với mức ý nghĩa alpha = 0.0001 các biến giải thích có ý nghĩa thống kê đối với mô hình bao gồm: age, overall, potential, international_reputation, positionGoalkeeper, crossing, finishing, short_passing, volleys, stamina, release_clause. Điều này cho thấy các biến này đều có ý nghĩa thống kê rất cao
- Việc có nhiều biến có ý nghĩa thống kê cao thường cho thấy mô hình hồi quy được xây dựng là phù hợp và có khả năng giải thích tốt dữ liệu. Mặc dù các biến này đều có ý nghĩa thống kê, cần cân nhắc về ý nghĩa thực tiễn của chúng. Một số biến có thể không có ý nghĩa thực tế mạnh dù có ý nghĩa thống kê cao.
1.2 Đánh giá mô hình
set.seed(21)
ctrl <- trainControl(method = "cv", number = 10)
cv_value_model <- train(value ~ ., data = simplified_data, method = "lm", trControl = ctrl)
print(cv_value_model)## Linear Regression
##
## 16643 samples
## 52 predictor
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 14978, 14980, 14977, 14978, 14980, 14979, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 0.5730172 0.9901131 0.2563575
##
## Tuning parameter 'intercept' was held constant at a value of TRUENhận xét:
- Dựa vào kết quả từ đoạn code trên, với các giá trị RRMSE = 0.575626, R-squared = 0.573021, và MAE = 0.2563513, so sánh với kết quả của hàm summary.lm trong mục Mô hình hồi quy tuyến tính cho cột value, ta nhận thấy mô hình có khả năng giải thích rất cao khi sử dụng tất cả các biến độc lập. Điều này được thể hiện qua giá trị R-squared gần bằng 1, cho thấy mô hình phù hợp với dữ liệu và có độ chính xác cao trong việc dự đoán.
- Tuy nhiên, mô hình hiện tại sử dụng rất nhiều biến giải thích, do đó, ta tìm cách đơn giản hóa mô hình
2. Lựa chọn mô hình
2.1 Hồi quy từng bước và cross-validation
predict.regsubsets <- function(object, newdata, id_model){
form <- as.formula(object$call[[2]])
x_mat <- model.matrix(form, newdata)
coef_est <- coef(object, id = id_model)
x_vars <- names(coef_est)
res <- x_mat[, x_vars] %*% coef_est
return(as.numeric(res))
}
# Chia bộ dữ liệu thành các folds
nrows <- nrow(simplified_data)
k <- 10
set.seed(21)
folds <- sample(rep(1:k,length=nrows))
cv_error_rj <- matrix(0,nrow=k,ncol=ncol(simplified_data))
for(r in 1:k){
data_train_r <- simplified_data[folds != r,]
data_test_r <- simplified_data[folds == r,]
steps_wise <- regsubsets(x=value ~ . ,data=data_train_r, method="forward", really.big = TRUE, nvmax = ncol(simplified_data))
for(j in 1:ncol(simplified_data)){
pred_rj <- predict(steps_wise, newdata = data_test_r,id_model = j)
cv_error_rj[r,j]<-sqrt(mean((data_test_r$value - pred_rj)^2))
}
}
cv_error <- colMeans(cv_error_rj)
ggplot(data= data.frame(x= c(1:ncol(simplified_data)),y=cv_error), mapping=aes(x=x,y=y)) +
geom_point() +
geom_line()+
labs(x="Number of predictors",y="RMSE") +
theme_bw()
Nhận xét:
- Ta lựa chọn số biến cho mô hình tối ưu có số biến là 12
steps_wise <- regsubsets(x = value~ ., data=simplified_data, method="forward", really.big = TRUE, nvmax = ncol(simplified_data))
selected_steps_wise <- coef(steps_wise, id = 12)
selected_steps_wise## (Intercept) age overall
## -0.630166142 -0.018947474 0.031994567
## potential wage international_reputation
## -0.019646451 0.006282471 0.401685586
## crossing finishing volleys
## -0.002041508 -0.003374028 0.003115535
## fk_accuracy stamina sliding_tackle
## 0.001716851 0.004330145 -0.002808977
## release_clause
## 0.485488417# Tạo công thức cho mô hình hồi quy
steps_wise_formula <- as.formula(paste("value ~", paste(names(coef(steps_wise, id = 12))[-1], collapse = " + ")))
steps_wise_formula <- update(steps_wise_formula, . ~ . - positionForward - positionMidfielder)
steps_wise_formula <- update(steps_wise_formula, . ~ . + position)
steps_wise_formula## value ~ age + overall + potential + wage + international_reputation +
## crossing + finishing + volleys + fk_accuracy + stamina +
## sliding_tackle + release_clause + position# Sử dụng k-folds cross-validation để đánh giá mô hình
set.seed(21)
ctrl <- trainControl(method = "cv", number = 10)
cv_steps_wise_value_model <- train(steps_wise_formula, data = simplified_data, method = "lm", trControl = ctrl)
print(cv_steps_wise_value_model)## Linear Regression
##
## 16643 samples
## 13 predictor
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 14978, 14980, 14977, 14978, 14980, 14979, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 0.5736017 0.9900932 0.2543194
##
## Tuning parameter 'intercept' was held constant at a value of TRUENhận xét:
- Sau khi giảm số biến giải thích xuống còn 12 biến, mô hình vẫn có khả năng giải thích rất cao với MSE = 0.5736017 và MAE = 0.2543194. Điều này cho thấy mô hình vẫn giữ được độ chính xác cao trong việc dự đoán giá trị của biến phụ thuộc, mặc dù số lượng biến giải thích đã giảm đáng kể. Việc giảm số biến giải thích giúp đơn giản hóa mô hình, làm cho nó dễ hiểu hơn và giảm nguy cơ overfitting (quá khớp) khi áp dụng trên dữ liệu mới.
2.2 Hồi quy lasso
Lasso (Least Absolute Shrinkage and selection Operator) nổi tiếng với khả năng lựa chọn biến bằng cách co một số hệ số về 0. Điều này rất hữu ích khi xử lý dữ liệu có nhiều biến dự đoán. Tuy nhiên, Lasso vẫn giả định mối quan hệ tuyến tính giữa các biến dự đoán còn lại và biến kết quả.
x_data <- model.matrix(value ~ ., data = simplified_data)[,-1]
y_data <- simplified_data$value
# Tìm hệ só lambda tối ưu bằng phương pháp cross-validation
set.seed(24)
out_cv_lasso_value <- cv.glmnet(x = x_data, y = y_data, alpha = 1, type.measure = "mse", nfolds = 10, family = "gaussian")
lambda_lasso_value <- out_cv_lasso_value$lambda.min
cat("Hệ số lambda tối ưu là: ", lambda_lasso_value)## Hệ số lambda tối ưu là: 0.02579507# Từ giá trị lambda tối ưu tìm được, fit mô hình trên toàn bộ tập data
out_lasso_md <- glmnet(x = x_data, y = y_data, alpha = 1, lambda = lambda_lasso_value, family = "gaussian")
# Các hệ số trong mô hình
coeff_lasso_md <- predict(out_lasso_md, type = "coefficients")
nonzero_coeff_lasso_md <- setNames(coeff_lasso_md@x, rownames(coeff_lasso_md)[coeff_lasso_md@i + 1])
nonzero_coeff_lasso_md## (Intercept) overall wage
## -1.1648293821 0.0134804901 0.0055903711
## international_reputation volleys stamina
## 0.3351279874 0.0006692361 0.0004724311
## release_clause
## 0.4863211056Nhận xét
- Các biến dự đoán có ý nghĩa: Các hệ số của các biến dự đoán sau không bằng không, cho thấy chúng có ý nghĩa trong việc dự đoán biến phản hồi:
overall: 0.0131,wage: 5.6743,international_reputation: 0.3363,volleys: 0.0007,reactions: 0.0002,stamina: 0.0005,release_clause: 0.4862- Diễn giải:
wagecó hệ số dương lớn nhất, cho thấy nó có tác động mạnh mẽ đến biến phản hồi.release_clausecũng có hệ số dương đáng kể, cho thấy nó là một biến dự đoán quan trọng.- Các hệ số không bằng không của
overall,international_reputation,volleys,reactions, vàstaminachỉ ra rằng các biến này cũng đóng góp vào mô hình, mặc dù ít hơn.
Đánh giá mô hình
set.seed(21)
k <- 10
folds <- sample(rep(1:k, length.out = nrow(simplified_data)))
cv_rmse <- numeric(k)
cv_rsquare <- numeric(k)
cv_rsquare_adjusted <- numeric(k)
for (i in 1:k) {
train_indices <- which(folds != i)
test_indices <- which(folds == i)
train_data <- simplified_data[train_indices, ]
test_data <- simplified_data[test_indices, ]
x_train <- model.matrix(value ~ ., data = train_data)[, -1]
y_train <- train_data$value
x_test <- model.matrix(value ~ ., data = test_data)[, -1]
y_test <- test_data$value
out_lasso <- glmnet(x = x_train, y = y_train, alpha = 1, lambda = lambda_lasso_value, family = "gaussian")
pred <- predict(out_lasso, s = lambda_lasso_value, newx = x_test)
cv_rmse[i] <- rmse(y_test, pred)
cv_rsquare[i] <- cor(y_test, pred)^2
cv_rsquare_adjusted[i] <- 1 - (1 - cv_rsquare[i]) * (length(y_test) - 1) / (length(y_test) - length(out_lasso$df))
}
avg_rmse <- mean(cv_rmse)
avg_rsquare <- mean(cv_rsquare)
avg_rsquare_adjusted <- mean(cv_rsquare_adjusted)
cat("Average RMSE: ", avg_rmse)## Average RMSE: 0.5781515cat("\nAverage R-squared: ", avg_rsquare)##
## Average R-squared: 0.9897874cat("\nAverage R-squared adjusted: ", avg_rsquare_adjusted)##
## Average R-squared adjusted: 0.98978743 Chuẩn đoán mô hình
3.1 Kiểm tra tính tuyến tính và Tính đồng nhất phương sai
pred_lasso_md <- predict(out_lasso_md, newx = x_data)
residual_lasso_md <- pred_lasso_md - simplified_data$value
ggplot(data = simplified_data, mapping = aes(x = pred_lasso_md, y = residual_lasso_md)) +
geom_point() +
geom_smooth(method = "loess", se = FALSE) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(x = "Fittted values", y = "Residuals") +
theme_minimal()## `geom_smooth()` using formula = 'y ~ x'
Nhận xét
- Đường trung bình phần dư (đường màu xanh) khá tương đồng với đường thẳng. Điều này có thể là dấu hiệu của mối quan hệ tuyến tính giữa value và một hoặc nhiều biến hồi quy.
- Các cầu thủ có giá trị thấp chiếm số lượng lớn, gây ra sự tập trung dày đặc trong đồ thị phần dư ở khu vực fitted values nhỏ. Ngược lại, số lượng cầu thủ có giá trị cao ít hơn, dẫn đến phần dư phân tán nhiều hơn ở các giá trị fitted values lớn. Điều này cho thấy phương sai của phần dư không đồng đều, vi phạm giả định về homoscedasticity (phương sai phần dư đồng đều) trong hồi quy tuyến tính.
3.2 Kiểm tra tính tuyến tính từng phần
fitted_lasso_md <- predict(out_lasso_md, newx = x_data)
resid_lasso_md <- y_data - fitted_lasso_md
selected_var_lasso_md <- names(nonzero_coeff_lasso_md[-1, drop = FALSE])
for(col_name in selected_var_lasso_md){
terms <- x_data[, col_name] * nonzero_coeff_lasso_md[col_name]
p <- ggplot(x_data, mapping = aes(x_data[, col_name], resid_lasso_md + terms)) +
geom_point() +
labs(x = col_name, y = "Partial Residuals") +
geom_smooth(method = "loess", se = FALSE, linetype = "dashed", color = "forestgreen") +
geom_line(aes(x = x_data[, col_name], y = terms), color = "blue")
theme_bw()
print(p)
}## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
Nhận xét
overall:
- Mối quan hệ phi tuyến tính: Đường smooth có xu hướng tăng lên khi giá trị “overall” tăng lên. Điều này cho thấy mối quan hệ giữa biến “overall” và biến kết quả có thể không tuyến tính - Heteroscedasticity (phương sai không đồng nhất): Sự phân tán của các điểm quanh đường smooth không đồng đều. Khi giá trị “overall” tăng lên (khoảng 75 trở lên), các điểm trở nên phân tán rộng hơn, đặc biệt là ở các giá trị trên 80.international_reputation:
- Mối quan hệ tuyến tính yếu: “international_reputation” trong mô hình Lasso. Đường màu xanh dương có độ dốc rất nhỏ, gần như bằng phẳng, cho thấy mối quan hệ tuyến tính giữa biến này và biến kết quả là rất yếu.
- Heteroscedasticity: Sự phân tán của các điểm có vẻ không đồng đều giữa các mức độ. Có vẻ như ở mức độ 3, các điểm phân tán rộng hơn so với các mức độ khác.
Voleys:
- Mối quan hệ phi tuyến tính nhẹ: Đường smooth cho thấy một mối quan hệ phi tuyến tính nhẹ giữa “volleys” và biến kết quả. Ở các giá trị “volleys” thấp và trung bình, đường smooth khá bằng phẳng, nhưng có xu hướng hơi tăng ở các giá trị “volleys” cao. - Phân bố dữ liệu tập trung: Phần lớn dữ liệu tập trung ở các giá trị “volleys” từ khoảng 0 đến 80. release_clause:
- Mối quan hệ tuyến tính mạnh mẽ: Đường smooth gần như trùng khớp với đường thẳng màu xanh dương, cho thấy một mối quan hệ tuyến tính rất mạnh mẽ giữa “release_clause” và biến kết quả. - Sự phân tán tương đối đồng đều: Các điểm dữ liệu phân tán khá đều quanh đường thẳng, không có dấu hiệu rõ ràng của heteroscedasticity (phương sai không đồng nhất).
4. Mở rộng mô hình
4.1 Hồi quy đa thức
- Dựa vào các biểu đồ kiểm tra tính tuyến tính từng phần, ta thấy rằng
các biến
overall,wage,volleysvàrelease_clausecó mối quan hệ phi tuyến tính với biến phụ thuộc. Do đó, ta sẽ xây dựng mô hình hồi quy đa thức cho các biến này.
poly_value_md <- lm(value ~ release_clause + poly(overall, 3) + poly(wage, 3) + international_reputation + poly(volleys,2) + stamina + release_clause*wage, data = simplified_data)
summary(poly_value_md)##
## Call:
## lm(formula = value ~ release_clause + poly(overall, 3) + poly(wage,
## 3) + international_reputation + poly(volleys, 2) + stamina +
## release_clause * wage, data = simplified_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.7527 -0.1113 -0.0088 0.1034 10.8366
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.335e-01 3.072e-02 7.600 3.12e-14 ***
## release_clause 4.051e-01 1.509e-03 268.366 < 2e-16 ***
## poly(overall, 3)1 8.089e+01 1.462e+00 55.331 < 2e-16 ***
## poly(overall, 3)2 7.567e+01 1.323e+00 57.214 < 2e-16 ***
## poly(overall, 3)3 3.804e+01 8.405e-01 45.255 < 2e-16 ***
## poly(wage, 3)1 -5.363e+01 1.909e+00 -28.101 < 2e-16 ***
## poly(wage, 3)2 -5.573e+01 1.532e+00 -36.378 < 2e-16 ***
## poly(wage, 3)3 -6.201e+00 5.747e-01 -10.790 < 2e-16 ***
## international_reputation -1.051e-01 1.626e-02 -6.463 1.06e-10 ***
## poly(volleys, 2)1 6.634e+00 6.413e-01 10.345 < 2e-16 ***
## poly(volleys, 2)2 1.264e+01 6.449e-01 19.602 < 2e-16 ***
## stamina 5.924e-03 3.522e-04 16.821 < 2e-16 ***
## wage NA NA NA NA
## release_clause:wage 3.704e-04 9.624e-06 38.483 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5115 on 16630 degrees of freedom
## Multiple R-squared: 0.992, Adjusted R-squared: 0.992
## F-statistic: 1.721e+05 on 12 and 16630 DF, p-value: < 2.2e-16Đánh giá mô hình bằng phương pháp cross-validation
set.seed(21)
k <- 10
folds <- sample(rep(1:k, length.out = nrow(simplified_data)))
cv_rmse <- numeric(k)
cv_rsquare <- numeric(k)
cv_rsquare_adjusted <- numeric(k)
for (i in 1:k) {
train_indices <- which(folds != i)
test_indices <- which(folds == i)
train_data <- simplified_data[train_indices, ]
test_data <- simplified_data[test_indices, ]
y_test <- test_data$value
poly_value_md <- lm(value ~ release_clause + poly(overall, 3) + poly(wage, 3) + international_reputation + poly(volleys,2) + stamina , data = train_data)
pred <- predict(poly_value_md, newdata = test_data)
cv_rmse[i] <- rmse(y_test, pred)
cv_rsquare[i] <- cor(y_test, pred)^2
cv_rsquare_adjusted[i] <- 1 - (1 - cv_rsquare[i]) * (length(y_test) - 1) / (length(y_test) - length(out_lasso$df))
}
avg_rmse <- mean(cv_rmse)
avg_rsquare <- mean(cv_rsquare)
avg_rsquare_adjusted <- mean(cv_rsquare_adjusted)
cat("Average RMSE: ", avg_rmse)## Average RMSE: 0.5401948cat("\nAverage R-squared: ", avg_rsquare)##
## Average R-squared: 0.9911411cat("\nAverage R-squared adjusted: ", avg_rsquare_adjusted)##
## Average R-squared adjusted: 0.99114114.2 Mô hình hồi quy tổng quát
gam_value_md <- gam(value ~ release_clause + s(overall) + s(wage) + international_reputation + s(volleys) + stamina + release_clause*wage, data = simplified_data)
summary(gam_value_md)##
## Family: gaussian
## Link function: identity
##
## Formula:
## value ~ release_clause + s(overall) + s(wage) + international_reputation +
## s(volleys) + stamina + release_clause * wage
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.014e+00 5.044e-02 20.101 < 2e-16 ***
## release_clause 4.083e-01 1.563e-03 261.254 < 2e-16 ***
## international_reputation -1.302e-01 1.603e-02 -8.124 4.83e-16 ***
## stamina 4.692e-03 3.679e-04 12.754 < 2e-16 ***
## wage -7.031e-02 4.268e-03 -16.471 < 2e-16 ***
## release_clause:wage 3.197e-04 1.134e-05 28.196 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(overall) 8.932 8.997 474.8 <2e-16 ***
## s(wage) 7.932 8.140 223.5 <2e-16 ***
## s(volleys) 8.978 9.000 115.6 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Rank: 32/33
## R-sq.(adj) = 0.993 Deviance explained = 99.3%
## GCV = 0.241 Scale est. = 0.24054 n = 16643predictions <- predict(gam_value_md, newdata = simplified_data)
mse <- mean((simplified_data$value - predictions)^2)
rmse <- sqrt(mse)
print(rmse)## [1] 0.4899834Đánh giá mô hình bằng phương pháp cross-validation
set.seed(21)
k <- 10
folds <- sample(rep(1:k, length.out = nrow(simplified_data)))
cv_rmse <- numeric(k)
cv_rsquare <- numeric(k)
cv_rsquare_adjusted <- numeric(k)
for (i in 1:k) {
train_indices <- which(folds != i)
test_indices <- which(folds == i)
train_data <- simplified_data[train_indices, ]
test_data <- simplified_data[test_indices, ]
y_test <- test_data$value
gam_value_md <- gam(value ~ release_clause + s(overall) + s(wage) + international_reputation + s(volleys) + stamina + release_clause*wage, data = train_data)
pred <- predict(gam_value_md, newdata = test_data)
cv_rmse[i] <- rmse(y_test, pred)
cv_rsquare[i] <- cor(y_test, pred)^2
cv_rsquare_adjusted[i] <- 1 - (1 - cv_rsquare[i]) * (length(y_test) - 1) / (length(y_test) - length(out_lasso$df))
}
avg_rmse <- mean(cv_rmse)
avg_rsquare <- mean(cv_rsquare)
avg_rsquare_adjusted <- mean(cv_rsquare_adjusted)
cat("Average RMSE: ", avg_rmse)## Average RMSE: 0.5166711cat("\nAverage R-squared: ", avg_rsquare)##
## Average R-squared: 0.9917792cat("\nAverage R-squared adjusted: ", avg_rsquare_adjusted)##
## Average R-squared adjusted: 0.99177925. Kết luận
- Sau khi xây dựng mô hình hồi quy tuyến tính, hồi quy lasso, hồi quy đa thức và hồi quy GAM, ta nhận thấy rằng mô hồi quy GAM có khả năng giải thích tốt nhất với Average RMSE: 0.5166711 Average R-squared: 0.9917792, Average R-squared adjusted: 0.9917792.
- Tổng số biến giải thích sử dụng cho mô hình này là 6 biến giải
thích, bao gồm:
release_clause,overall,wage,international_reputation,volleysvàstamina.
IV. Linear Regresion cho wage
wage_lm <- lm(wage ~ ., data = simplified_data)
summary(wage_lm)##
## Call:
## lm(formula = wage ~ ., data = simplified_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -154.819 -2.072 -0.239 1.581 195.650
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.183e+02 1.566e+02 4.587 4.53e-06 ***
## age 2.732e-01 4.095e-02 6.671 2.62e-11 ***
## overall -3.915e-02 4.690e-02 -0.835 0.403933
## potential 2.165e-02 3.402e-02 0.637 0.524399
## value 2.232e+00 1.443e-01 15.471 < 2e-16 ***
## preferred_footRight 3.728e-02 2.086e-01 0.179 0.858152
## international_reputation 9.678e+00 2.995e-01 32.317 < 2e-16 ***
## weak_foot -1.681e-01 1.379e-01 -1.219 0.222699
## skill_moves -9.454e-01 2.256e-01 -4.190 2.80e-05 ***
## attacking_work_rateLow -4.702e-01 4.538e-01 -1.036 0.300182
## attacking_work_rateMedium 5.003e-03 2.178e-01 0.023 0.981674
## deffensive_work_rate Low -1.197e+00 3.941e-01 -3.037 0.002394 **
## deffensive_work_rate Medium -9.628e-01 2.449e-01 -3.932 8.47e-05 ***
## body_typeNormal -3.246e-01 1.942e-01 -1.671 0.094706 .
## body_typeStocky -8.245e-01 3.971e-01 -2.076 0.037880 *
## positionForward 1.254e+00 4.772e-01 2.627 0.008614 **
## positionGoalkeeper 4.857e+00 1.799e+00 2.700 0.006937 **
## positionMidfielder -2.817e-01 3.358e-01 -0.839 0.401551
## jersey_number 1.412e-02 5.465e-03 2.585 0.009758 **
## joined -2.529e-01 4.155e-02 -6.085 1.19e-09 ***
## contract_valid_until -1.167e-01 6.944e-02 -1.680 0.092934 .
## height 1.357e-01 6.489e-02 2.091 0.036537 *
## weight 1.939e-02 9.822e-03 1.974 0.048385 *
## crossing 3.842e-02 1.186e-02 3.241 0.001194 **
## finishing -2.069e-02 1.439e-02 -1.437 0.150666
## heading_accuracy 1.683e-02 1.267e-02 1.328 0.184249
## short_passing -1.246e-02 2.030e-02 -0.614 0.539470
## volleys -9.718e-03 1.229e-02 -0.791 0.429220
## dribbling 1.899e-02 1.778e-02 1.068 0.285683
## curve 4.183e-03 1.213e-02 0.345 0.730257
## fk_accuracy -3.701e-02 1.077e-02 -3.436 0.000591 ***
## long_passing -2.351e-02 1.494e-02 -1.573 0.115647
## ball_control 4.969e-02 2.216e-02 2.242 0.024945 *
## acceleration -1.537e-02 1.664e-02 -0.924 0.355628
## sprint_speed 2.035e-02 1.544e-02 1.318 0.187380
## agility -1.337e-05 1.235e-02 -0.001 0.999136
## reactions -3.647e-02 1.883e-02 -1.937 0.052771 .
## balance 3.492e-02 1.262e-02 2.766 0.005678 **
## shot_power 2.591e-02 1.252e-02 2.069 0.038567 *
## jumping 1.110e-03 8.728e-03 0.127 0.898806
## stamina -4.060e-02 1.042e-02 -3.897 9.78e-05 ***
## strength -2.356e-02 1.227e-02 -1.920 0.054919 .
## long_shots 1.825e-03 1.327e-02 0.137 0.890667
## aggression -6.621e-04 9.482e-03 -0.070 0.944329
## interceptions 1.345e-02 1.352e-02 0.995 0.319690
## positioning 2.545e-03 1.328e-02 0.192 0.848007
## vision 4.666e-03 1.240e-02 0.376 0.706754
## penalties 3.875e-02 1.167e-02 3.321 0.000899 ***
## composure -1.655e-02 1.382e-02 -1.198 0.231026
## marking -1.717e-02 1.101e-02 -1.560 0.118850
## standing_tackle -8.809e-03 2.015e-02 -0.437 0.661933
## sliding_tackle 7.356e-02 1.876e-02 3.920 8.88e-05 ***
## gk_diving 5.206e-04 2.591e-02 0.020 0.983971
## gk_handling 3.398e-02 2.609e-02 1.303 0.192762
## gk_kicking -1.186e-02 2.415e-02 -0.491 0.623374
## gk_positioning -5.471e-02 2.518e-02 -2.173 0.029793 *
## gk_reflexes -6.646e-03 2.561e-02 -0.260 0.795239
## release_clause 3.560e-01 7.228e-02 4.925 8.51e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.67 on 16585 degrees of freedom
## Multiple R-squared: 0.7711, Adjusted R-squared: 0.7703
## F-statistic: 980.4 on 57 and 16585 DF, p-value: < 2.2e-16- Ta nhận thấy độ chính xác của mô hình rất thấp
1. Lựa chọn mô hình
1.1 Hồi quy từng bước và cross-validation
predict.regsubsets <- function(object, newdata, id_model){
form <- as.formula(object$call[[2]])
x_mat <- model.matrix(form, newdata)
coef_est <- coef(object, id = id_model)
x_vars <- names(coef_est)
res <- x_mat[, x_vars] %*% coef_est
return(as.numeric(res))
}
# Chia bộ dữ liệu thành các folds
nrows <- nrow(simplified_data)
k <- 10
set.seed(23)
folds <- sample(rep(1:k,length=nrows))
cv_error_rj <- matrix(0,nrow=k,ncol=ncol(simplified_data))
for(r in 1:k){
data_train_r <- simplified_data[folds != r,]
data_test_r <- simplified_data[folds == r,]
steps_wise <- regsubsets(x=wage ~ . ,data=data_train_r, method="forward", really.big = TRUE, nvmax = ncol(simplified_data))
for(j in 1:ncol(simplified_data)){
pred_rj <- predict(steps_wise, newdata = data_test_r,id_model = j)
cv_error_rj[r,j]<-sqrt(mean((data_test_r$value - pred_rj)^2))
}
}
cv_error <- colMeans(cv_error_rj)
ggplot(data= data.frame(x= c(1:ncol(simplified_data)),y=cv_error), mapping = aes(x = x,y=y)) +
geom_point() +
geom_line()+
labs(x="Number of predictors",y="RMSE") +
theme_bw()
out_subset_wage <- regsubsets(x = wage ~ ., data = simplified_data, method = "forward", nvmax = ncol(simplified_data))
coef(out_subset_wage, id = which.min(cv_error))## (Intercept) value
## 1.428039 3.352893Đánh giá mô hình
# Tạo công thức cho mô hình hồi quy
subset_wage_formula <- as.formula(paste("wage ~", paste(names(coef(out_subset_wage, id = which.min(cv_error)))[-1], collapse = " + ")))
#subset_wage_formula <- update(subset_wage_formula, . ~ . - positionMidfielder)
#subset_wage_formula <- update(subset_wage_formula, . ~ . + position)
subset_wage_formula## wage ~ value# Sử dụng k-folds cross-validation để đánh giá mô hình
set.seed(21)
ctrl <- trainControl(method = "cv", number = 10)
cv_subset_wage <- train(subset_wage_formula, data = simplified_data, method = "lm", trControl = ctrl)
print(cv_subset_wage)## Linear Regression
##
## 16643 samples
## 1 predictor
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 14978, 14980, 14978, 14979, 14980, 14979, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 11.27782 0.736994 4.767674
##
## Tuning parameter 'intercept' was held constant at a value of TRUE1.2 Hồi quy lasso
Lasso (Least Absolute Shrinkage and selection Operator) nổi tiếng với khả năng lựa chọn biến bằng cách co một số hệ số về 0. Điều này rất hữu ích khi xử lý dữ liệu có nhiều biến dự đoán. Tuy nhiên, Lasso vẫn giả định mối quan hệ tuyến tính giữa các biến dự đoán còn lại và biến kết quả.
x_data_wage <- model.matrix(wage ~ ., data = simplified_data)[,-1]
y_data_wage <- simplified_data$wage
# Tìm hệ só lambda tối ưu bằng phương pháp cross-validation
set.seed(21)
out_cv_lasso_wage <- cv.glmnet(x = x_data_wage, y = y_data_wage, alpha = 1, type.measure = "mse", nfolds = 10, family = "gaussian")
lambda_lasso_wage <- out_cv_lasso_wage$lambda.min
cat("Hệ số lambda tối ưu là: ", lambda_lasso_wage)## Hệ số lambda tối ưu là: 0.01353117# Từ giá trị lambda tối ưu tìm được, fit mô hình trên toàn bộ tập data
out_lasso_wage <- glmnet(x = x_data_wage, y = y_data_wage, alpha = 1, lambda = lambda_lasso_wage, family = "gaussian")
# Các hệ số trong mô hình
coeff_lasso_wage <- predict(out_lasso_wage, type = "coefficients")
nonzero_coeff_lasso_wage <- setNames(coeff_lasso_wage@x, rownames(coeff_lasso_wage)[coeff_lasso_wage@i + 1])
sort(nonzero_coeff_lasso_wage, decreasing = T)## (Intercept) international_reputation
## 6.958199e+02 9.720065e+00
## positionGoalkeeper value
## 2.720660e+00 2.321381e+00
## positionForward release_clause
## 1.081515e+00 3.056566e-01
## age height
## 2.357464e-01 1.136092e-01
## sliding_tackle crossing
## 6.288031e-02 3.665628e-02
## ball_control penalties
## 3.308576e-02 3.252654e-02
## balance shot_power
## 2.718679e-02 2.111946e-02
## weight dribbling
## 1.539182e-02 1.525324e-02
## gk_handling jersey_number
## 1.392103e-02 1.384133e-02
## heading_accuracy interceptions
## 1.033243e-02 5.719712e-03
## sprint_speed curve
## 3.291576e-03 1.283051e-04
## jumping acceleration
## 9.939096e-05 -3.373238e-04
## short_passing volleys
## -7.539697e-04 -3.246271e-03
## finishing marking
## -1.013665e-02 -1.106648e-02
## composure overall
## -1.118500e-02 -1.766929e-02
## strength long_passing
## -1.879400e-02 -1.953088e-02
## gk_positioning reactions
## -2.442599e-02 -2.965283e-02
## fk_accuracy stamina
## -3.175789e-02 -3.775269e-02
## contract_valid_until weak_foot
## -1.083770e-01 -1.463178e-01
## joined body_typeNormal
## -2.483883e-01 -2.715708e-01
## positionMidfielder attacking_work_rateLow
## -4.079299e-01 -4.415127e-01
## body_typeStocky deffensive_work_rate Medium
## -7.319435e-01 -8.779270e-01
## skill_moves deffensive_work_rate Low
## -8.874190e-01 -1.075323e+00Nhận xét
- Từ hệ số của các biến trong mô hình hồi quy Lasso, ta thấy chỉ có hệ số của các biến giải thích như
overall,international_reputation,volleys,reactions,stamina,release_clausekhác 0, cho thấy các biến này có ảnh hưởng đến biến kết quảwage.- Từ hệ số của các biến trong mô hình hồi quy Lasso, ta các biến giải thích sau international_reputation, position, value, release_clause có hệ số lớn nhât
Đánh giá mô hình
set.seed(21)
k <- 10
folds <- sample(rep(1:k, length.out = nrow(simplified_data)))
cv_rmse <- numeric(k)
cv_rsquare <- numeric(k)
cv_rsquare_adjusted <- numeric(k)
for (i in 1:k) {
train_indices <- which(folds != i)
test_indices <- which(folds == i)
train_data <- simplified_data[train_indices, ]
test_data <- simplified_data[test_indices, ]
x_train <- model.matrix(wage ~ ., data = train_data)[, -1]
y_train <- train_data$wage
x_test <- model.matrix(wage ~ ., data = test_data)[, -1]
y_test <- test_data$wage
out_lasso <- glmnet(x = x_train, y = y_train, alpha = 1, lambda = lambda_lasso_wage, family = "gaussian")
pred <- predict(out_lasso, s = lambda_lasso_wage, newx = x_test)
cv_rmse[i] <- rmse(y_test, pred)
cv_rsquare[i] <- cor(y_test, pred)^2
cv_rsquare_adjusted[i] <- 1 - (1 - cv_rsquare[i]) * (length(y_test) - 1) / (length(y_test) - length(out_lasso$df))
}
avg_rmse <- mean(cv_rmse)
avg_rsquare <- mean(cv_rsquare)
avg_rsquare_adjusted <- mean(cv_rsquare_adjusted)
cat("Average RMSE: ", avg_rmse)## Average RMSE: 10.664cat("\nAverage R-squared: ", avg_rsquare)##
## Average R-squared: 0.7649655cat("\nAverage R-squared adjusted: ", avg_rsquare_adjusted)##
## Average R-squared adjusted: 0.7649655- Từ kết quả của mô hình lasso, ta đánh giá mô hình chỉ dựa trên các biến giải thích có hệ số lớn như: international_reputation, position, value, deffensive_work_rate, skill_moves
set.seed(21)
k <- 10
folds <- sample(rep(1:k, length.out = nrow(simplified_data)))
cv_rmse <- numeric(k)
cv_rsquare <- numeric(k)
cv_rsquare_adjusted <- numeric(k)
for (i in 1:k) {
train_indices <- which(folds != i)
test_indices <- which(folds == i)
train_data <- simplified_data[train_indices, ]
test_data <- simplified_data[test_indices, ]
x_train <- model.matrix(wage ~ international_reputation + position + value + deffensive_work_rate + skill_moves, data = train_data)[, -1]
y_train <- train_data$wage
x_test <- model.matrix(wage ~ international_reputation + position + value + deffensive_work_rate + skill_moves, data = test_data)[, -1]
y_test <- test_data$wage
out_lasso <- glmnet(x = x_train, y = y_train, alpha = 1, lambda = lambda_lasso_wage, family = "gaussian")
pred <- predict(out_lasso, s = lambda_lasso_wage, newx = x_test)
cv_rmse[i] <- rmse(y_test, pred)
cv_rsquare[i] <- cor(y_test, pred)^2
cv_rsquare_adjusted[i] <- 1 - (1 - cv_rsquare[i]) * (length(y_test) - 1) / (length(y_test) - length(out_lasso$df))
}
avg_rmse <- mean(cv_rmse)
avg_rsquare <- mean(cv_rsquare)
avg_rsquare_adjusted <- mean(cv_rsquare_adjusted)
cat("Average RMSE: ", avg_rmse)## Average RMSE: 10.74997cat("\nAverage R-squared: ", avg_rsquare)##
## Average R-squared: 0.7610392cat("\nAverage R-squared adjusted: ", avg_rsquare_adjusted)##
## Average R-squared adjusted: 0.7610392Nhận xét
- Mô hình Lasso sử dụng các biến giải thích với hệ số khác 0 cho kết quả: Average RMSE: 10.66493, Average R-squared: 0.7649227, Average R-squared adjusted: 0.7649227. Điều này cho thấy mô hình không có khả năng giải thích tốt
- Sau khi chọn ra các biến giải thích như: international_reputation, position, value, deffensive_work_rate, skill_moves, mô hình có khả năng giải thích tốt hơn với Average RMSE: 10.74997, Average R-squared: 0.7610392, Average R-squared adjusted: 0.7610392. Ta thấy cả 2 mô hình không có sự khác biệt lớn về độ chính xác giữa chúng. Do đó, ta sẽ chọn mô hình Lasso sử dụng ít biến giải thích hơn để giảm độ phức tạp của mô hình.
x_data_wage <- model.matrix(wage ~ international_reputation + position + value + deffensive_work_rate + skill_moves, data = simplified_data)[,-1]
y_data_wage <- simplified_data$wage
out_cv_lasso_wage <- cv.glmnet(x = x_data_wage, y = y_data_wage, alpha = 1, type.measure = "mse", nfolds = 10, family = "gaussian")
lambda_lasso_wage <- out_cv_lasso_wage$lambda.min
out_lasso_wage <- glmnet(x = x_data_wage, y = y_data_wage, alpha = 1, lambda = lambda_lasso_wage, family = "gaussian")
coeff_lasso_wage <- predict(out_lasso_wage, type = "coefficients")
nonzero_coeff_lasso_wage <- setNames(coeff_lasso_wage@x, rownames(coeff_lasso_wage)[coeff_lasso_wage@i + 1])
nonzero_coeff_lasso_wage## (Intercept) international_reputation
## -6.9870834 10.8537193
## positionForward positionGoalkeeper
## -0.4267046 -1.1183289
## positionMidfielder value
## -1.1458339 2.8598020
## deffensive_work_rate Low deffensive_work_rate Medium
## -1.6097663 -1.2082737
## skill_moves
## -0.34568512. Chuẩn đoán mô hình
2.1 Kiểm tra tính tuyến tính và Tính đồng nhất phương sai
pred_lasso_wage <- predict(out_lasso_wage, newx = x_data_wage)
residual_lasso_wage <- pred_lasso_wage - simplified_data$wage
ggplot(data = simplified_data, mapping = aes(x = pred_lasso_wage, y = residual_lasso_wage)) +
geom_point() +
geom_smooth(method = "loess", se = FALSE) +
geom_hline(yintercept = 0, linetype = "dashed") +
labs(x = "Fittted values", y = "Residuals") +
theme_minimal()## `geom_smooth()` using formula = 'y ~ x'
Nhận xét
- Đường trung bình phần dư (đường màu xanh) có xu hướng cong xuống khi giá trị biến phụ thuộc tăng lên. Điều này có thể là dấu hiệu của mối quan hệ phi tuyến tuyến tính giữa wage và các biến giải thích
- Các cầu thủ có giá trị thấp chiếm số lượng lớn, gây ra sự tập trung dày đặc trong đồ thị phần dư ở khu vực fitted values nhỏ. Ngược lại, số lượng cầu thủ có giá trị cao ít hơn, dẫn đến phần dư phân tán nhiều hơn ở các giá trị fitted values lớn. Ở các cầu thủ có giá trị lớn, các điểm phân tán xa đường ngang residuals 0 hơn. Vì vậy, mô hình không đáp ứng được giả định về tính đồng nhất phương sai.
2.2 Kiểm tra tính tuyến tính từng phần
fitted_lasso_wage <- predict(out_lasso_wage, newx = x_data_wage)
resid_lasso_wage <- simplified_data$wage - fitted_lasso_wage
for(col_name in colnames(x_data_wage)){
terms <- x_data_wage[, col_name] * nonzero_coeff_lasso_wage[col_name]
p <- ggplot(x_data_wage, mapping = aes(x_data_wage[, col_name], resid_lasso_wage + terms)) +
geom_point() +
labs(x = col_name, y = "Partial Residuals") +
geom_smooth(method = "loess", se = FALSE, linetype = "dashed", color = "forestgreen") +
geom_line(aes(x = x_data_wage[, col_name], y = terms), color = "blue")
theme_bw()
print(p)
}## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
Nhận xét
international_reputation:
- Mối quan hệ tuyến tính mạnh mẽ: Đường smooth gần như trùng khớp với đường thẳng màu xanh dương, cho thấy một mối quan hệ tuyến tính rất mạnh mẽ giữa “international_reputation” và biến kết quả.
positionForward, positionGoalkeeper, positionMidfieler, defensive_work_rate, skill_moves:
- Các phần dư có sự phân tán lớn, đặc biệt ở cả hai đầu. Điều này có thể cho thấy một mức độ biến thiên đáng kể mà biến positionForward, positionaMildfielder không thể giải thích.
3. Mở rộng mô hình
3.1 Mô hình hồi quy đa thức
poly_wage_md <- lm(wage ~ international_reputation + poly(value, 3) + poly(skill_moves, 3) , data = simplified_data)
summary(poly_wage_md)##
## Call:
## lm(formula = wage ~ international_reputation + poly(value, 3) +
## poly(skill_moves, 3), data = simplified_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -163.952 -1.906 -0.955 0.844 191.303
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.4788 0.3234 -10.756 < 2e-16 ***
## international_reputation 11.7479 0.2804 41.897 < 2e-16 ***
## poly(value, 3)1 2122.5095 15.2559 139.128 < 2e-16 ***
## poly(value, 3)2 145.6990 11.1851 13.026 < 2e-16 ***
## poly(value, 3)3 -35.8573 10.9297 -3.281 0.00104 **
## poly(skill_moves, 3)1 -19.5235 11.8291 -1.650 0.09887 .
## poly(skill_moves, 3)2 -107.2838 11.3589 -9.445 < 2e-16 ***
## poly(skill_moves, 3)3 -83.2415 10.8345 -7.683 1.64e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 10.71 on 16635 degrees of freedom
## Multiple R-squared: 0.7689, Adjusted R-squared: 0.7688
## F-statistic: 7905 on 7 and 16635 DF, p-value: < 2.2e-163.2 Mô hình tuyến tính tổng quát
gam_wage_md <- gam(wage ~ international_reputation + s(value) + value * international_reputation + skill_moves, data = simplified_data)
summary(gam_wage_md)##
## Family: gaussian
## Link function: identity
##
## Formula:
## wage ~ international_reputation + s(value) + value * international_reputation +
## skill_moves
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.09120 1.10488 12.754 < 2e-16 ***
## international_reputation 3.35370 0.36048 9.303 < 2e-16 ***
## value -4.59308 0.41214 -11.145 < 2e-16 ***
## skill_moves -0.35949 0.11642 -3.088 0.00202 **
## international_reputation:value 0.91223 0.02533 36.014 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(value) 8.159 8.175 164.2 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Rank: 13/14
## R-sq.(adj) = 0.789 Deviance explained = 79%
## GCV = 104.46 Scale est. = 104.38 n = 16643predictions <- predict(gam_wage_md, newdata = simplified_data)
mse <- mean((simplified_data$wage - predictions)^2)
rmse <- sqrt(mse)
print(rmse)## [1] 10.212754. Kết luận
- So với mô hình ban đầu, mô hình cuối cùng chỉ sử dụng các biến international_reputation, value, skill_moves có kết quả tương đồng, nhưng mô hình cuối cùng đơn giản hơn rất nhiều
- Kết quả từ mô hình với R-squared-adjusted = 0.789, và RMSE = 10.21275 cho thấy mô hình hồi quy cho biến giải thích wage không có khả năng giải thích cho bộ dữ liệu